{$L- CURVE_FIT.LIST}

{
                                                  program by:  EDWARD J GRUNDLER
      Changes made by reviewer:
          Reversed ranges in type SCREEN.
          Replaced PAGE(OUTPUT) with CLEARSCREEN.
          Replaced definitions of RLF and CEOL with those in SCREEN_STUFF.
          Replaced missing unit EDS_STUFF with REAL_INPUT.
                   ___Henry E. Baumgarten
}

PROGRAM CURVE_FIT;
USES {$U GRAPH.LIBRARY} SCREEN_STUFF, REAL_INPUT;

CONST MAXIMUM_POINTS = 50;
      MAXIMUM_FUNCTIONS = 10;
      MAX_PLUS_1 = 11;
TYPE MATRIX = ARRAY[1..MAXIMUM_FUNCTIONS,1..MAX_PLUS_1] OF REAL;
     REAL_ARRAY = ARRAY[1..MAXIMUM_POINTS] OF REAL;
     SCREEN = PACKED ARRAY[0..239,0..319] OF BOOLEAN;
VAR ERROR_SUM,SQUARE_SUM:REAL;
    C,X,Y,YLS,ERROR,ERROR_SQ:REAL_ARRAY;
    R,NUMB_OF_POINTS,NUMBER_OF_FUNCTIONS:INTEGER;
    TITLE:STRING;
    ANSWER,TYPE_FUNCTION,NEW:CHAR;
    ERROR_FLAG,STARTED:BOOLEAN;
    A:MATRIX;

PROCEDURE DRAWLINE(VAR R:INTEGER; VAR S:SCREEN; RW,XS,YS,DX,DY,INK:INTEGER);
EXTERNAL;

PROCEDURE QUERY(Q:STRING);
BEGIN
  WRITELN(Q);
  REPEAT READ(KEYBOARD,ANSWER) UNTIL ANSWER IN ['Y','y','N','n']
END;

FUNCTION F(X:REAL;I:INTEGER):REAL;
  FUNCTION F1:REAL;
  BEGIN
    I:=I-1;
    IF I=0 THEN F1:=1 ELSE F1:=F1*X
  END;
BEGIN {F}
  CASE TYPE_FUNCTION OF
  'P','p':F:=F1;
  END {OF CASE}
END;

PROCEDURE DISPLAY;
VAR I:INTEGER;
    XL,XU,YL,YU,XRNG,YRNG,XINC,YINC:REAL;
    S:SCREEN;

  PROCEDURE SET_SCREEN;
  VAR J,XSTART,YSTART:INTEGER;
  BEGIN {SET_SCREEN}
    CLEARSCREEN;
    FILLCHAR(S,SIZEOF(S),CHR(0));
    UNITWRITE(3,S,63);
    FOR J:=0 TO 9 DO DRAWLINE(R,S,20,24,44+J*20,5,0,1);
    FOR J:=0 TO 10 DO DRAWLINE(R,S,20,39+J*28,229,0,-5,1);
    GOTOXY(0,4);
    FOR J:=0 TO 9 DO
      BEGIN
        WRITELN(YU-J*YRNG/9:2:2);
        IF J<>9 THEN WRITELN
      END;
    FOR J:=0 TO 10 DO
      BEGIN
        GOTOXY(3+7*J,23);
        WRITE(XL+J*XRNG/10:7:2)
      END;
    GOTOXY(0,0);
    IF ABS(319-XU/XINC)>319 THEN XSTART:=0 ELSE XSTART:=ROUND(319-XU/XINC);
    IF XSTART IN [24..319] THEN DRAWLINE(R,S,20,XSTART,39,0,190,1);
    IF ABS(44+YU/YINC)>239 THEN YSTART:=0 ELSE YSTART:=ROUND(44+YU/YINC);
    IF YSTART IN [40..229] THEN DRAWLINE(R,S,20,24,YSTART,295,0,1)
  END; {SET_SCREEN}
  
  PROCEDURE PLOT_POINTS;
  VAR XSTART,YSTART,J:INTEGER;
  BEGIN {PLOT_POINTS}
    FOR J:=1 TO NUMB_OF_POINTS DO
      BEGIN
        XSTART:=319-ROUND((XU-X[J])/XINC);
        YSTART:=ROUND((YU-Y[J])/YINC)+44;
        DRAWLINE(R,S,20,XSTART-5,YSTART+5,5,-5,1);
        DRAWLINE(R,S,20,XSTART-5,YSTART-5,5,5,1)
      END
  END; {PLOT_POINTS}
  
  PROCEDURE DRAW_CURVE;
  VAR J,K,YSTART,YSTOP:INTEGER;
      X1,Y1,X2,Y2:REAL;
  BEGIN {DRAW_CURVE}
    FOR J:=39 TO 318 DO
      BEGIN
        X1:=XL+(J-39)*XINC;
        X2:=XL+(J-38)*XINC;
        Y1:=0;
        FOR K:=1 TO NUMBER_OF_FUNCTIONS DO Y1:=Y1+C[K]*F(X1,K);
        Y2:=0;
        FOR K:=1 TO NUMBER_OF_FUNCTIONS DO Y2:=Y2+C[K]*F(X2,K);
        IF ABS(44+(YU-Y1)/YINC)>239 THEN YSTART:=0
        ELSE YSTART:=ROUND(44+(YU-Y1)/YINC);
        IF ABS(44+(YU-Y2)/YINC)>239 THEN YSTOP:=0
        ELSE YSTOP:=ROUND(44+(YU-Y2)/YINC);
        IF (YSTART IN [40..229]) AND (YSTOP IN [40..229]) THEN
        DRAWLINE(R,S,20,J,YSTART,1,YSTOP-YSTART,1)
      END
  END; {DRAW_CURVE}

BEGIN {DISPLAY}
  XL:=X[1];
  XU:=X[1];
  YL:=Y[1];
  YU:=Y[1];
  FOR I:=2 TO NUMB_OF_POINTS DO
    BEGIN
      IF X[I]<XL THEN XL:=X[I] ELSE IF X[I]>XU THEN XU:=X[I];
      IF Y[I]<YL THEN YL:=Y[I] ELSE IF Y[I]>YU THEN YU:=Y[I]
    END;
  IF (XU-XL)=0 THEN
    BEGIN
      XL:=XL-0.05;
      XU:=XU+0.05
    END
  ELSE
    BEGIN
      XL:=XL-(XU-XL)/18;
      XU:=XU+(XU-XL)/19
    END;
  XRNG:=XU-XL;
  IF (YU-YL)=0 THEN
    BEGIN
      YU:=YU+0.045;
      YL:=YL-0.045
    END
  ELSE
    BEGIN
      YL:=YL-(YU-YL)/18;
      YU:=YU+(YU-YL)/19
    END;
  YRNG:=YU-YL;
  XINC:=XRNG/280;
  YINC:=YRNG/180;
  SET_SCREEN;
  PLOT_POINTS;
  DRAW_CURVE;
  WRITE('press RETURN to continue');
  READLN;
  UNITWRITE(3,S,7)
END; {DISPLAY}
      
PROCEDURE OUTPUT_LEAST_SQUARES;
VAR I:INTEGER;
BEGIN
  FOR I:=1 TO NUMB_OF_POINTS DO
    BEGIN
      IF (I MOD 20)=1 THEN
        BEGIN
          CLEARSCREEN;
          WRITELN(' N       X(N)           Y(N)      Y=L.S.APPROX.      ',
                  'ERR=Y-Y(N)      ERR^2');
          WRITELN
        END;
      WRITELN(I:2,X[I]:15:6,Y[I]:15:6,YLS[I]:15:6,ERROR[I]:15:6,
              ERROR_SQ[I]:15:6);
      IF ((I MOD 20)=0) OR (I=NUMB_OF_POINTS) THEN
        BEGIN
          WRITE('press RETURN to continue');
          IF I=NUMB_OF_POINTS THEN
            WRITE('                  SUMS=',ERROR_SUM:15:6,SQUARE_SUM:15:6);
          READLN
        END
    END
END;

PROCEDURE LEAST_SQUARES;
VAR I,J:INTEGER;
BEGIN
  ERROR_SUM:=0;
  SQUARE_SUM:=0;
  FOR I:=1 TO NUMB_OF_POINTS DO
    BEGIN
      WRITE('.');
      YLS[I]:=0;
      FOR J:=1 TO NUMBER_OF_FUNCTIONS DO
        YLS[I]:=YLS[I]+C[J]*F(X[I],J);
      ERROR[I]:=YLS[I]-Y[I];
      ERROR_SQ[I]:=SQR(ERROR[I]);
      ERROR_SUM:=ERROR_SUM+ERROR[I];
      SQUARE_SUM:=SQUARE_SUM+ERROR_SQ[I]
    END
END;

PROCEDURE OUTPUT_C;
VAR I:INTEGER;
BEGIN
  CLEARSCREEN;
  WRITELN(TITLE);
  WRITELN;
  FOR I:=1 TO NUMBER_OF_FUNCTIONS DO
    WRITELN('C(',I:2,') = ',C[I]:12:6);
  WRITELN;
  WRITELN('press RETURN to continue')
END;

PROCEDURE GAUSS_SEIDEL;
VAR I,J,K:INTEGER;
    X,ERROR:REAL;
BEGIN {GAUSS_SEIDEL}
  K:=0;
  REPEAT
    ERROR:=0;
    K:=K+1;
    WRITE('.');
    FOR I:=1 TO NUMBER_OF_FUNCTIONS DO
      BEGIN
        X:=A[I,NUMBER_OF_FUNCTIONS+1];
        FOR J:=1 TO NUMBER_OF_FUNCTIONS DO IF NOT (J=I) THEN X:=X-A[I,J]*C[J];
        X:=X/A[I,I];
        ERROR:=ERROR+ABS(X-C[I]);
        C[I]:=X
      END;
    IF (K MOD 50)=0 THEN WRITELN
  UNTIL (ERROR=0) OR (K=150);
  FOR I:=1 TO NUMBER_OF_FUNCTIONS DO A[I,NUMBER_OF_FUNCTIONS+1]:=C[I]
END; {GAUSS_SEIDEL}

PROCEDURE ROW_ECHELON(A:MATRIX; ROWS,COLS:INTEGER);
TYPE INTEGER_ARRAY = ARRAY[1..MAX_PLUS_1] OF INTEGER;
VAR I,J,K:INTEGER;
    SCALAR:INTEGER_ARRAY;
  PROCEDURE CLEAR_COL(I,J,N:INTEGER);
  VAR K,L,M:INTEGER;
      MULT:REAL;
  BEGIN {CLEAR_COL}
    FOR K:=I TO J DO
      BEGIN
        IF A[K,N]<>0 THEN
          BEGIN
            MULT:=-A[K,N]/A[N,N];
            FOR L:=1 TO COLS DO
              A[K,L]:=A[K,L]+MULT*A[N,L]
          END
      END
  END; {CLEAR_COL}
  PROCEDURE SWAP_ROWS(I:INTEGER);
  VAR J,L,M:INTEGER;
      HOLD,MAX:REAL;
  BEGIN {SWAP_ROWS}
    MAX:=ABS(A[I,I]);
    M:=I;
    FOR J:=I+1 TO ROWS DO
      BEGIN
        IF ABS(A[J,I])>MAX THEN
          BEGIN
            MAX:=ABS(A[J,I]);
            M:=J
          END
      END;
    IF M<>I THEN
      BEGIN
        FOR L:=1 TO COLS DO
          BEGIN
            HOLD:=A[I,L];
            A[I,L]:=A[M,L];
            A[M,L]:=HOLD
          END
      END
  END; {SWAP_ROWS}
  PROCEDURE DIVIDE_ROW(I:INTEGER);
  VAR M:INTEGER;
      DIVISOR:REAL;
  BEGIN {DIVIDE_ROW}
    DIVISOR:=A[I,I];
    FOR M:=1 TO COLS DO A[I,M]:=A[I,M]/DIVISOR
  END; {DIVIDE_ROW}
  PROCEDURE UNSCALE;
  VAR I,J:INTEGER;
      X:REAL;
  BEGIN {UNSCALE}
    FOR J:=1 TO NUMBER_OF_FUNCTIONS+1 DO
      BEGIN
        X:=1;
        FOR I:=1 TO SCALAR[J] DO X:=X*2;
        FOR I:=1 TO NUMBER_OF_FUNCTIONS DO A[I,J]:=A[I,J]*X
      END
  END; {UNSCALE}
  PROCEDURE SCALE;
  VAR I,J:INTEGER;
      X:REAL;
  BEGIN {SCALE}
    FOR I:=1 TO NUMBER_OF_FUNCTIONS+1 DO SCALAR[I]:=0;
    FOR J:=1 TO NUMBER_OF_FUNCTIONS+1 DO
      BEGIN
        X:=ABS(A[1,J]);
        FOR I:=2 TO NUMBER_OF_FUNCTIONS DO IF ABS(A[I,J])>X THEN X:=ABS(A[I,J]);
        WHILE X>10 DO
          BEGIN
            SCALAR[J]:=SCALAR[J]+1;
            X:=X/2
          END;
        X:=1;
        FOR I:=1 TO SCALAR[J] DO X:=X*2;
        FOR I:=1 TO NUMBER_OF_FUNCTIONS DO A[I,J]:=A[I,J]/X
      END
  END; {SCALE}
BEGIN {ROW_ECHELON}
  SCALE;
  FOR I:=1 TO ROWS DO
    BEGIN
      WRITE('.');
      SWAP_ROWS(I);
      IF A[I,I]<>0 THEN CLEAR_COL(I+1,ROWS,I)
      ELSE
        BEGIN
          ERROR_FLAG:=TRUE;
          EXIT(ROW_ECHELON)
        END
    END;
  UNSCALE;
  FOR I:=ROWS DOWNTO 1 DO
    BEGIN
      WRITE('.');
      DIVIDE_ROW(I);
      CLEAR_COL(1,I-1,I)
    END;
  FOR I:=1 TO ROWS DO C[I]:=A[I,NUMBER_OF_FUNCTIONS+1];
  WRITELN
END; {ROW_ECHELON}

PROCEDURE B_SOLVE;
VAR I,J:INTEGER;
BEGIN
  FOR I:=1 TO NUMBER_OF_FUNCTIONS DO
    BEGIN
      WRITE('.');
      A[I,NUMBER_OF_FUNCTIONS+1]:=0;
      FOR J:=1 TO NUMB_OF_POINTS DO
        A[I,NUMBER_OF_FUNCTIONS+1]:=A[I,NUMBER_OF_FUNCTIONS+1]+F(X[J],I)*Y[J]
    END;
  WRITELN
END;

PROCEDURE A_SOLVE;
VAR I,J,K:INTEGER;
BEGIN
  FOR I:=1 TO NUMBER_OF_FUNCTIONS DO
    BEGIN
      FOR J:=1 TO NUMBER_OF_FUNCTIONS DO
        BEGIN
          WRITE('.');
          A[I,J]:=0;
          FOR K:=1 TO NUMB_OF_POINTS DO A[I,J]:=A[I,J]+F(X[K],I)*F(X[K],J)
        END;
      WRITELN
    END
END;

PROCEDURE GET_HOW_MANY;
BEGIN
  REPEAT
    CLEARSCREEN;
    WRITELN('ENTER THE NUMBER OF FUNCTIONS');
    WRITELN;
    WRITELN('FOR EXAMPLE;');
    WRITELN('IF YOU ARE SOLVING FOR A POLYNOMIAL');
    WRITELN('ENTER 3 FOR C( 1)+C( 2)X+C( 3)X^2');
    REPEAT;
      WRITE('NUMBER OF FUNCTIONS = ');
      NUMBER_OF_FUNCTIONS:=ROUND(INPUT_VALUE);
      IF NUMBER_OF_FUNCTIONS>MAXIMUM_FUNCTIONS THEN
        WRITE('MAXIMUM IS CURRENTLY ',MAXIMUM_FUNCTIONS,' ')
      ELSE
        IF NUMBER_OF_FUNCTIONS>NUMB_OF_POINTS THEN
          WRITE('ALL FUNCTIONS GREATER THAN ',NUMB_OF_POINTS,
                ' WOULD EVALUATE TO 0.  ');
    UNTIL NUMBER_OF_FUNCTIONS IN[1..NUMB_OF_POINTS];
    WRITE('NUMBER OF FUNCTIONS IS ',NUMBER_OF_FUNCTIONS);
    QUERY(' -- IS THAT CORRECT')
  UNTIL ANSWER IN ['Y','y']
END;

PROCEDURE VERIFY_DATA;
VAR I:INTEGER;
  PROCEDURE FIX_DATA;
  VAR X_OR_Y:CHAR;
      J:INTEGER;
  BEGIN {FIX_DATA}
    REPEAT
      REPEAT
        WRITELN(RLF,'DO YOU WANT TO CHANGE AN X OR Y VALUE? ',CEOL);
        READ(KEYBOARD,X_OR_Y)
      UNTIL X_OR_Y IN ['X','x','Y','y'];
      REPEAT
        WRITE(RLF,'WHICH ',X_OR_Y,' DO YOU WANT TO CHANGE? ',CEOL);
        J:=ROUND(INPUT_VALUE)
      UNTIL J IN [((I-1) DIV 20*20+1)..I];
      WRITE(RLF,'ENTER VALUE FOR ',X_OR_Y,'(',J,') ',CEOL);
      IF X_OR_Y IN ['X','x'] THEN
        BEGIN
          X[J]:=INPUT_VALUE;
          GOTOXY(2,(J-1) MOD 20+2);
          WRITE(X[J]:15:3)
        END
      ELSE
        BEGIN
          Y[J]:=INPUT_VALUE;
          GOTOXY(17,(J-1) MOD 20+2);
          WRITE(Y[J]:10:3)
        END;
      GOTOXY(0,(I-1) MOD 20+3);
      QUERY('ARE THERE ANY MORE POINTS TO CHANGE?')
    UNTIL ANSWER IN ['N','n']
  END; {FIX_DATA}
BEGIN {VERIFY_DATA}
  FOR I:=1 TO NUMB_OF_POINTS DO
    BEGIN
      IF (I MOD 20)=1 THEN
        BEGIN
          CLEARSCREEN;
          WRITELN(' I          X         Y');
          WRITELN;
        END;
      WRITELN(I:2,X[I]:15:3,Y[I]:10:3);
      IF (I=NUMB_OF_POINTS) OR ((I MOD 20)=0) THEN
        BEGIN
          QUERY('ARE ALL OF THE ABOVE DATA CORRECT?');
          IF ANSWER IN ['N','n'] THEN FIX_DATA
        END
    END
END; {VERIFY_DATA}

PROCEDURE GET_DATA;
VAR I:INTEGER;
BEGIN
  CLEARSCREEN;
  REPEAT
    REPEAT
      WRITE('HOW MANY DATA POINTS ARE THERE? ');
      NUMB_OF_POINTS:=ROUND(INPUT_VALUE);
      IF NUMB_OF_POINTS>MAXIMUM_POINTS THEN
        WRITELN('CURRENT MAXIMUM IS ',MAXIMUM_POINTS,' POINTS')
    UNTIL NUMB_OF_POINTS IN [1..MAXIMUM_POINTS];
    WRITE('NUMBER OF POINTS IS ',NUMB_OF_POINTS);
    QUERY(' -- IS THAT CORRECT?')
  UNTIL ANSWER IN ['Y','y'];
  FOR I:=1 TO NUMB_OF_POINTS DO
    BEGIN
      WRITE('ENTER THE VALUE FOR X(',I,') ');
      X[I]:=INPUT_VALUE;
      WRITE('ENTER THE VALUE FOR Y(',I,') ');
      Y[I]:=INPUT_VALUE
    END
END;

PROCEDURE NEW_OR_OLD;
BEGIN
  IF STARTED THEN
    BEGIN
      CLEARSCREEN;
      WRITE('do you want to enter N(ew data or use the O(ld data? ');
      REPEAT READ(KEYBOARD,NEW) UNTIL NEW IN ['N','n','O','o']
    END
  ELSE
    BEGIN
      STARTED:=TRUE;
      NEW:='N'
    END
END;

PROCEDURE GET_TYPE;
BEGIN
  CLEARSCREEN;
  WRITELN('THE TYPE FUNCTIONS AVAILABLE ARE:');
  WRITELN;
  WRITELN('P(olynomial');
  WRITELN;
  WRITELN('PLEASE SELECT');
  REPEAT READ(KEYBOARD,TYPE_FUNCTION) UNTIL TYPE_FUNCTION IN ['P','p']
END;

PROCEDURE GET_TITLE;
BEGIN
  REPEAT
    WRITELN('ENTER THE TITLE OF THIS CURVE ');
    READLN(TITLE);
    WRITELN(TITLE);
    QUERY('IS THAT CORRECT?')
  UNTIL ANSWER IN ['Y','y']
END;

PROCEDURE INITIALIZE;
BEGIN
  DEFINESTRINGS;
  CLEARSCREEN;
  STARTED:=FALSE;
  WRITELN('THIS PROGRAM ACCEPTS DATA POINTS, P(X,Y) AS INPUT');
  WRITELN('AND GIVES A LEAST SQUARES APPROXIMATION CURVE');
  WRITELN('FOR THE DATA.');
  WRITELN;
  WRITELN('SEVERAL OPTIONS WILL BE PRESENTED AT');
  WRITELN('APPROPRIATE TIMES.  WHEN THESE OPTIONS');
  WRITELN('ARE PRESENTED, ENTER ONLY THE UPPER CASE');
  WRITELN('LETTER.  FOR EXAMPLE:  ENTER P FOR');
  WRITELN('P(olynomial.');
  WRITELN
END;

BEGIN
  INITIALIZE;
  REPEAT
    ERROR_FLAG:=FALSE;
    GET_TITLE;
    GET_TYPE;
    NEW_OR_OLD;
    IF NEW IN ['N','n'] THEN
      BEGIN
        GET_DATA;
        VERIFY_DATA
      END;
    GET_HOW_MANY;
    A_SOLVE;
    B_SOLVE;
    ROW_ECHELON(A,NUMBER_OF_FUNCTIONS,NUMBER_OF_FUNCTIONS+1);
    GAUSS_SEIDEL;
    IF ERROR_FLAG THEN
      BEGIN
        WRITELN('THERE IS NO UNIQUE SOLUTION');
        WRITELN;
        WRITELN('press RETURN to continue');
        READLN
      END
    ELSE
      BEGIN
        OUTPUT_C;
        LEAST_SQUARES;
        READLN;
        OUTPUT_LEAST_SQUARES
      END;
    DISPLAY;
    CLEARSCREEN;
    QUERY('IS THAT THE LAST CURVE FIT?')
  UNTIL ANSWER IN ['Y','y']
END.
